setwd("~/Desktop/redistributing power/CPS submission")
land <- read.csv("taiwan_land.csv")

install.packages("ggplot2")
library(ggplot2)
############### variable transformation #############

#ethnic fractionalization
land$quan_per <- land$quanzhou_jpcensus / land$pop_jpcensus
land$zhang_per <- land$zhangzhou_jpcensus / land$pop_jpcensus
land$hakka_per <- land$guangdong_jpcensus / land$pop_jpcensus
land$otherhan_per <- land$otherhan_census / land$pop_jpcensus
land$settlement_frac <- 1 - (land$quan_per^2) - (land$zhang_per^2) - (land$hakka_per^2) - (land$otherhan_per^2)

#tenancy
land$tenancy_jp <- (land$tenant_jp + land$mixed_jp) / (land$independent_jp + land$mixed_jp + land$tenant_jp)
land$independency_jp <- (land$independent_jp) / (land$independent_jp + land$mixed_jp + land$tenant_jp)
land$tenancy_57 <- (land$tenant_57) / (land$tenant_57 + land$mixed_57 + land$independent_57)

#membership purges
land$sponsor_ratio <- land$sponsor_members_56 / (land$postreform_member + land$missing_member_56)
land$disqualify_ratio <- land$disqualified_member_53 / land$prereform_member_53
############## farming association cooptation
######grassroot officers
#1956
a <- land$member_rep_re_56 / (land$member_rep_re_56 + land$member_rep_new_56)
b <- land$team_cap_re_56 / (land$team_cap_re_56 + land$team_cap_new_56)
land$officer_reelect_56 <- (a + b) / 2
#1959
a.2 <- land$member_rep_re_59 / (land$member_rep_re_59 + land$member_rep_new_59)
b.2 <- land$team_cap_re_59 / (land$team_cap_re_59 + land$team_cap_new_59)
land$officer_reelect_59 <- (a.2 + b.2) / 2
#1962
a.3 <- land$member_rep_re_62 / (land$member_rep_re_62 + land$member_rep_new_62)
b.3 <- land$team_cap_re_62 / (land$team_cap_re_62 + land$team_cap_new_62)
land$officer_reelect_62 <- (a.3 + b.3) / 2
#1965
a.4 <- land$member_rep_re_65 / (land$member_rep_re_65 + land$member_rep_new_65)
b.4 <- land$team_cap_re_65 / (land$team_cap_re_65 + land$team_cap_new_65)
land$officer_reelect_65 <- (a.4 + b.4) / 2
#aggregate
land$officer_reelect <- (rowSums(land[,c("officer_reelect_56", "officer_reelect_59", "officer_reelect_62", "officer_reelect_65")], na.rm=TRUE))/4
mean(land$member_rep_re_56 + land$member_rep_new_56 + land$member_rep_re_59 + land$member_rep_new_59 + land$member_rep_re_62 + land$member_rep_new_62 + land$member_rep_re_65 + land$member_rep_new_65 , na.rm = T)

######board
#1956
c <- land$councilor_re_56 / (land$councilor_re_56 + land$councilor_new_56)
d <- land$auditor_re_56 / (land$auditor_re_56 + land$auditor_new_56)
land$board_reelect_56 <- (c + d) / 2
#1959
c.2 <- land$councilor_re_59 / (land$councilor_re_59 + land$councilor_new_59)
d.2 <- land$auditor_re_59 / (land$auditor_re_59 + land$auditor_new_59)
land$board_reelect_59 <- (c.2 + d.2) / 2
#1962
c.3 <- land$councilor_re_62 / (land$councilor_re_62 + land$councilor_new_62)
d.3 <- land$auditor_re_62 / (land$auditor_re_62 + land$auditor_new_62)
land$board_reelect_62  <- (c.3 + d.3) / 2
#1965
c.4 <- land$councilor_re_65 / (land$councilor_re_65 + land$councilor_new_65)
d.4 <- land$auditor_re_65 / (land$auditor_re_65 + land$auditor_new_65)
land$board_reelect_65  <- (c.4 + d.4) / 2

#aggregate
land$board_reelect <- (rowSums(land[,c("board_reelect_56", "board_reelect_59", "board_reelect_62", "board_reelect_65")], na.rm=TRUE))/4

######exec
land$executive_reelect <- (land$exec_re_56 + land$exec_re_59 + land$exec_re_62 + land$exec_re_65 + land$chairman_re_56 + land$chairman_re_59 + land$chairman_re_62 + land$chairman_re_65 + land$auditorexec_re_56 + land$auditorexec_re_59 + land$auditorexec_re_62 + land$auditorexec_re_65) / 12

######## total tally
land$position_reelect <- (land$officer_reelect + land$board_reelect + land$executive_reelect) / 3

land$position_reelect_56 <- (land$exec_re_56 + land$board_reelect_56 + land$officer_reelect_56)/3
land$position_reelect_59 <- (land$exec_re_59 + land$board_reelect_59 + land$officer_reelect_59)/3
land$position_reelect_62 <- (land$exec_re_62 + land$board_reelect_62 + land$officer_reelect_62)/3
land$position_reelect_65 <- (land$exec_re_65 + land$board_reelect_65 + land$officer_reelect_65)/3

#average number of positions per township:
land$all_56 <- land$member_rep_re_56 + land$member_rep_new_56 + land$team_cap_re_56 + land$team_cap_new_56 + land$councilor_re_56 + land$councilor_new_56 + land$auditor_re_56 + land$auditor_new_56 + 1
land$all_59 <- land$member_rep_re_59 + land$member_rep_new_59 + land$team_cap_re_59 + land$team_cap_new_59 + land$councilor_re_59 + land$councilor_new_59 + land$auditor_re_59 + land$auditor_new_59 + 1
land$all_62 <- land$member_rep_re_62 + land$member_rep_new_62 + land$team_cap_re_62 + land$team_cap_new_62 + land$councilor_re_62 + land$councilor_new_62 + land$auditor_re_62 + land$auditor_new_62 + 1
land$all_65 <- land$member_rep_re_65 + land$member_rep_new_65 + land$team_cap_re_65 + land$team_cap_new_65 + land$councilor_re_65 + land$councilor_new_65 + land$auditor_re_65 + land$auditor_new_65 + 1

land$all <- land$all_56 + land$all_59 + land$all_62 + land$all_65 
mean(land$all, na.rm = T)

##### main models ##################################
#####################################################
################ revised submission August 9

####main results (Text Table 2)
m.1 <- lm(position_reelect ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + county_EN, data = land)
summary(m.1)

m.2 <- lm(sponsor_ratio ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + county_EN, data = land)
summary(m.2)

#####robustness check (Appendix Table A)
m.a.1 <- lm(officer_reelect ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + county_EN, data = land)
summary(m.a.1)

m.a.2 <- lm(position_reelect_56 ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + county_EN, data = land)
summary(m.a.2)

m.a.3 <- lm(position_reelect_59 ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + county_EN, data = land)
summary(m.a.3)

m.a.4 <- lm(position_reelect_62 ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + county_EN, data = land)
summary(m.a.4)

m.a.5 <- lm(position_reelect_65 ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + county_EN, data = land)
summary(m.a.5)

#####robustness check (Appendix Table B)
m.b.1 <- lm(disqualify_ratio ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_jp + county_EN, data = land)
summary(m.b.1)

m.b.2 <- lm(disqualify_ratio ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + county_EN, data = land)
summary(m.b.2)

######## predicted effects plot (Text Graph 2)
install.packages("ggeffects")
library(ggeffects)
m.1.p <- ggpredict(m.1, terms = "req_per")
plot_a <- plot(m.1.p) + xlab("Ratio of requisitioned land") + ylab("Ratio of 
reelected leaders") + ggtitle("")

m.2.p <- ggpredict(m.2, terms = "twhead")
plot_b <- plot(m.2.p) + xlab("Native representation under colonialism") + ylab(" Ratio of associate membership") + ggtitle("")

p_list <- list(plot_a, plot_b)
library(ggpubr)
figure <- ggarrange(plotlist=p_list,
                    ncol = 2, nrow = 1)


######## robustness check (Appendix Table C)
install.packages('betareg')
library(betareg)

m.c.1 <- betareg(position_reelect ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + county_EN, data = land, link = c("logit"))
summary(m.c.1)

m.c.2 <- betareg(sponsor_ratio ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + county_EN, data = land, link = c("logit"))
summary(m.c.2)

m.c.3 <- lm(position_reelect ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + prefecture_EN, data = land)
summary(m.c.3)

m.c.4 <- lm(sponsor_ratio ~ req_per + twhead + settlement_frac + colonizer_land_per + postreform_member + tenancy_57 + prefecture_EN, data = land)
summary(m.c.4)

stargazer(m.c.1, m.c.2, m.c.3, m.c.4, style = 'apsr', type = 'html')

######leadership continuity graph (Text Graph 3)
year <- c("1956", "1959", "1962", "1965")
board <- c(0.3782119, 0.4406845, 0.3297666, 0.3283003)
officer <- c(0.4748298, 0.5180235, 0.4500319, 0.4532945)
cycle <- data.frame(year, board, officer)

install.packages("reshape2") 
library(reshape2)
melt.cycle <- melt(cycle,id.vars="year")

ggplot(melt.cycle, aes(year, value, color = variable, group = variable)) + 
  geom_line() + xlab("Electoral cycle") + ylab("Reeelcted ratio") +
  ylim(0.2, 0.6) + guides(color = guide_legend(title = "Leadership level")) +
  scale_color_manual(values=c("red", "blue")) + theme(text = element_text(size = 15))


#########Geographical distribution of two main independent variables

tw_new <- readOGR(dsn = "~/Desktop/gis/TWN", layer = "Townships", encoding = "UTF-8")
tw_new.df <- fortify(tw_new, region = "UNI_ID")

library(dplyr)
land$id <- as.character(land$id)
final.plot <- left_join(tw_new.df, land, by = "id")

install.packages("mapproj")
library(mapproj)

map_1 <- ggplot() + 
  geom_polygon(data = final.plot, aes(x = long, y = lat, group = group, fill = position_reelect), size = 0.25) + 
  theme_void() +
  labs(title="Ratio of re-elected farmers association leaders, 1956-1965") +
  scale_fill_gradientn(colours = brewer.pal(9,"Blues")) 

map_2 <- ggplot() + 
  geom_polygon(data = final.plot, aes(x = long, y = lat, group = group, fill = sponsor_ratio), size = 0.25) + 
  theme_void() +
  labs(title="purged farmers association members, 1956") +
  scale_fill_gradientn(colours = brewer.pal(9,"Reds")) 

figure <- ggarrange(map_1, map_2)

  